This script generates static and animated maps of prescribed burns in Three Rivers parks. It aligns the master burn spreadsheet with the spatial burn unit layer before plotting the maps. The maps can be used in annual reports and public presentations.
# Loading required packages for analysis:
library(here)
library(tidyverse)
library(ggplot2)
library(sf)
library(ggspatial)
library(gganimate)
library(gifski)
library(transformr)
library(ggrepel)
source(here("Scripts", "1_Data_Import.R"))
source(here("Scripts", "2_Spatial_Data_Import.R"))
source(here("Scripts", "3_Data_Cleaning.R"))
source(here("Scripts", "4_Covariate_Data_Prep.R"))
# source(here("Scripts", "5_Prep_Species.R"))The burn units spatial data are imported with the script “2_Spatial_Data_Import.R”. The burn data are imported and processed in the script “4_Covariate_Data_Prep.R”. We make copies of both of those layers for further modification here. There is lots of manual work to align both.
# Load spatial data on burn units
units <- spatial.mgmt.burns
# Load burn dates
burn_dates <- burns %>%
select(Burn_unit_ID_new_2015, Burn_unit_description, year, burn_status_simple) %>%
mutate(burn_status_simple = ifelse(is.na(burn_status_simple), 0, burn_status_simple))
# Manually make burn unit IDs in burn_dates match unit IDs in units
burn_dates <- burn_dates %>%
mutate(Burn_unit_ID_new_2015 = ifelse(Burn_unit_description == "East Prairie", "NHTC 1", Burn_unit_ID_new_2015)) %>%
mutate(Burn_unit_ID_new_2015 = ifelse(Burn_unit_description == "West Prairie", "NHTC 2", Burn_unit_ID_new_2015)) %>%
mutate(Burn_unit_ID_new_2015 = ifelse(Burn_unit_ID_new_2015 == "REB 1", "REB 1A", Burn_unit_ID_new_2015)) %>%
mutate(Burn_unit_ID_new_2015 = ifelse(Burn_unit_ID_new_2015 == "REB 2", "REB 1B", Burn_unit_ID_new_2015)) %>%
mutate(Burn_unit_ID_new_2015 = ifelse(Burn_unit_description == "Open Area", "SPR NA", Burn_unit_ID_new_2015)) %>%
# Carver 22 has west and east (a & b) in the burn dates layer, but just one polygon in the spatial layer. CAR 22 B in the burn dates layer has no burns marked. so I just use the A burns and apply them to the whole 22 polygon (both halves).
filter(Burn_unit_ID_new_2015 != "CAR 22 B") %>% # drop 22 B
mutate(Burn_unit_ID_new_2015 = ifelse(Burn_unit_description == "Old field south of mitigation pond. West Half", "CAR 22", Burn_unit_ID_new_2015)) %>%
# Murphy 7 now has north and south in the burn dates layer (though both labeled 7B...), but just one polygon in the spatial layer. Going to combine both rows in the burns data and call both 7.
filter(Burn_unit_description != "1996 Metro Prairie-East - sw marsh. South Half") %>% # drop South half
mutate(burn_status_simple = ifelse(Burn_unit_ID_new_2015 == "MUH 7 B" & year == "2023", "1", burn_status_simple)) %>% # rescue the one burn in the south and apply to unit 7. TK this will need to be updated in future years (until spatial layer is updated...)
mutate(Burn_unit_ID_new_2015 = ifelse(Burn_unit_ID_new_2015 == "MUH 7 B", "MUH 7", Burn_unit_ID_new_2015))
# Join unit data (burn unit polygons) and burn data (data frame indicating when each unit was burned)
units <- units %>%
select(ABBREVIATI, BURN_ID_2015, NRM_UNIT, Unit_ID) %>%
rename(
park = ABBREVIATI,
id = BURN_ID_2015
) %>%
# Alter some unit names/numbers to resolve differences between burn units & burn dates
mutate(id = ifelse(park == "ELM_A", "3A", id)) %>%
mutate(id = ifelse(park == "ELM_B", "3B", id)) %>%
mutate(id = ifelse(park == "ELM_C", "3C", id)) %>%
mutate(park = ifelse(park %in% c("ELM_A", "ELM_B", "ELM_C"), "ELM", park)) %>%
mutate(id = ifelse(park == "GW" & NRM_UNIT == "GW1", "1", id)) %>%
mutate(id = ifelse(park == "GW" & NRM_UNIT == "GW2", "2", id)) %>%
mutate(id = ifelse(park == "GW" & NRM_UNIT == "GW3", "3", id)) %>%
mutate(id = ifelse(park == "HYL A" & id == "4", "4A", id)) %>%
mutate(id = ifelse(park == "HYL B" & id == "4", "4B", id)) %>%
mutate(id = ifelse(park == "HYL A" & id == "7", "7A", id)) %>%
mutate(id = ifelse(park == "HYL B" & id == "7", "7B", id)) %>%
mutate(id = ifelse(park == "HYL C" & id == "7", "7C", id)) %>%
mutate(id = ifelse(park == "HYL" & NRM_UNIT == "B" & id == "6", "6A", id)) %>% # mis-match in NRM unit letter and renamed value intentional, pending confirmation [TK]
mutate(id = ifelse(park == "HYL" & NRM_UNIT == "A" & id == "6", "6B", id)) %>% # mis-match in NRM unit letter and renamed value intentional, pending confirmation [TK]
mutate(id = ifelse(park == "HYL C" & id == "6", "6C", id)) %>%
mutate(park = ifelse(park %in% c("HYL A", "HYL B", "HYL C"), "HYL", park)) %>%
mutate(park = ifelse(park %in% c("MUL_E", "MUL_W"), "MUL", park)) %>%
mutate(id = ifelse(Unit_ID %in% c("NHTC1"), "1", id)) %>%
mutate(id = ifelse(Unit_ID %in% c("NHTC2"), "2", id)) %>%
mutate(id = ifelse(park == "REB A", "1A", id)) %>%
mutate(id = ifelse(park == "REB B", "1B", id)) %>%
mutate(park = ifelse(park %in% c("REB A", "REB B"), "REB", park)) %>%
mutate(unit = paste0(park, " ", id)) %>% # Make match how units listed in burn data
select(park, unit)
# Match units in spatial layer with names in date layer
units <- units %>%
left_join(burn_dates, by = c("unit" = "Burn_unit_ID_new_2015")) %>%
drop_na() %>%
arrange(unit, year)We use habitat types as a basemap over which to plot the burns.
# Generate basemap (habitat types)
hab_groups_forest <- c("Oak Woods", "Early Successional Forest", "Shrubland", "Maple Basswood", "Flood Plain Forest", "Tamarack Bog")
hab_groups_grasslands <- c("Oldfield", "Prairie")
hab_groups_wetlands <- c("Wetland Cattail", "Wetland Sedge", "Wetland Lily", "Beach Mudflat")
hab_group_water <- c("Open Water")
hab_groups_cultivated <- c("Planted Cultivated")
basemap_habs <- spatial.habitats.current %>%
mutate(type_simplified = ifelse(MLCCS_Data %in% hab_groups_forest, "forest",
ifelse(MLCCS_Data %in% hab_groups_grasslands, "grassland",
ifelse(MLCCS_Data %in% hab_groups_wetlands, "wetland",
ifelse(MLCCS_Data %in% hab_groups_cultivated, "cultivated",
ifelse(MLCCS_Data %in% hab_group_water, "water", NA)
)
)
)
)) %>%
filter(!is.na(type_simplified))
basemap_cols <- c("forest" = "darkseagreen", "grassland" = "lightgoldenrod", "cultivated" = "wheat3", "wetland" = "lightblue", "water" = "skyblue4")
basemap <- ggplot() +
geom_sf(data = basemap_habs, aes(fill = type_simplified), color = NA, show.legend = FALSE) +
scale_fill_manual(values = basemap_cols) +
theme_void()
basemap# Generate park name lookup table for parks with burns (relates codes in burn data to park names in park park boundaries dataset).
park_name_lookup <- data.frame(
park = c("BAK", "BRY", "CAR", "CRO", "ELM", "FR", "GW", "HYL", "KW", "MUH", "MUL", "NHTC", "REB", "SPR"),
Name = c("Baker", "Bryant", "Carver", "Crow-Hassan", "Elm Creek", "French", "Gale Woods", "Hyland", "Kingswood", "Murphy-Hanrehan", "The Landing", "NHTC", "Lake Rebecca", "Spring Lake Park")
)
# Function for quickly defining map extent in ggplot based on a particular park. Park name must match that listed in park boundaries sf.
map_extent_park <- function(park) {
row <- which(spatial.boundaries.parks$Name == park)
bbox <- st_bbox(spatial.boundaries.parks[row, ])
coord_sf(
xlim = c(bbox$xmin, bbox$xmax),
ylim = c(bbox$ymin, bbox$ymax)
)
}
# Function for mapping burns in a particular park/year
map_burns <- function(park, display_year) {
units_filter_year <- filter(units, year == display_year)
cols <- c("1" = "orangered3", "forest" = "darkseagreen", "grassland" = "lightgoldenrod", "wetland" = "lightblue", "water" = "skyblue4", cultivated = "wheat3")
ggplot() +
geom_sf(data = spatial.boundaries.parks, fill = "grey97", size = 0) +
geom_sf(data = basemap_habs, aes(fill = type_simplified), color = NA, show.legend = FALSE, alpha = 0.7) +
# geom_sf(data = spatial.infrastructure.trails %>%
# filter(st_geometry_type(.) %in% c("MULTILINESTRING")) %>%
# st_zm(drop = TRUE),
# linetype = "dotted", size = .5, alpha = .8, color = "tan4") +
geom_sf(
data = spatial.infrastructure.roads %>%
filter(st_geometry_type(.) %in% c("MULTILINESTRING")) %>%
st_zm(drop = TRUE),
linetype = "solid", size = .4, alpha = .8, color = "tan4"
) +
geom_sf(data = units_filter_year %>% filter(burn_status_simple == 1), aes(fill = as.factor(burn_status_simple)), show.legend = FALSE, color = "white", alpha = 0.75) +
geom_sf(data = spatial.boundaries.parks, fill = NA, size = 0.7) +
scale_fill_manual(values = cols) +
scale_y_continuous(breaks = seq(40, 50, by = 0.02)) +
scale_x_continuous(breaks = seq(-95, -90, by = 0.02)) +
theme_minimal() +
theme(
axis.title = element_blank(),
plot.title = element_text(size = 11, face = "bold", color = "orangered3"),
plot.subtitle = element_text(size = 10, face = "plain")
) +
map_extent_park(park) +
annotation_scale(location = "br", style = "bar", bar_cols = c("grey40", "grey40"), line_col = "grey40", text_col = "grey40", line_width = 0.5) +
labs(
title = "Prescribed burns",
subtitle = paste0(park, ": ", display_year)
)
}Example map, Crow-Hassan 1990:
maps_2023 <- lapply(park_name_lookup$Name, map_burns, "2023")
names(maps_2023) <- park_name_lookup$Name
maps_2023## $Baker
##
## $Bryant
##
## $Carver
##
## $`Crow-Hassan`
##
## $`Elm Creek`
##
## $French
##
## $`Gale Woods`
##
## $Hyland
##
## $Kingswood
##
## $`Murphy-Hanrehan`
##
## $`The Landing`
##
## $NHTC
##
## $`Lake Rebecca`
##
## $`Spring Lake Park`
# Function for generating an animation of burns in a particular park
map_burns_for_ani <- function(park) {
# Define colors
cols <- c("1" = "orangered3", "cultivated" = "wheat3", "forest" = "darkseagreen", "grassland" = "lightgoldenrod", "wetland" = "lightblue", "water" = "skyblue4")
# Pull burned untis only
units_burned <- units %>% filter(burn_status_simple == 1)
# Get park code from park name
park_code <- park_name_lookup$park[which(park_name_lookup$Name == park)]
# Calculate the minimum start year for the current park
min_start_year <- min(as.numeric(as.character(units_burned$year[which(units_burned$park == park_code)])))
# Create a subset of data starting from the year before the first burn
units_subset <- units_burned %>% filter(park == park_code, as.numeric(as.character(year)) >= min_start_year - 1)
# For mapping purposes, we need a fake/placeholder row in the birn units database (one for each unique combination of park and year). I'll have the geometry be a small circle centered on Minneapolis (out of frame of any park).
fake_geometry <- st_sfc(st_point(c(-93.2650, 44.9778)))
fake_geometry <- st_buffer(fake_geometry, dist = 100)
st_crs(fake_geometry) <- st_crs(units)
# Create a data frame with unique combinations of year and park
fake_data <- expand.grid(year = unique(units$year), park = unique(units$park), burn_status_simple = 1, unit = NA, Burn_unit_description = NA) %>%
filter(as.numeric(as.character(year)) >= min_start_year - 1) %>%
filter(park == park_code)
# Combine the data frame with the spatial data
fake_df <- st_sf(fake_data, Shape = fake_geometry)
# Bind the fake rows to the original units_burned dataframe
units_subset <- rbind(units_subset, fake_df)
# Generate basemap for animation
plot <- ggplot() +
geom_sf(data = spatial.boundaries.parks, fill = "grey97", size = 0) +
geom_sf(data = basemap_habs, aes(fill = type_simplified), color = NA, show.legend = FALSE, alpha = 0.7) +
# geom_sf(data = spatial.infrastructure.trails %>%
# filter(st_geometry_type(.) %in% c("MULTILINESTRING")) %>%
# st_zm(drop = TRUE),
# linetype = "dotted", size = .5, alpha = .8, color = "tan4") +
geom_sf(
data = spatial.infrastructure.roads %>%
filter(st_geometry_type(.) %in% c("MULTILINESTRING")) %>%
st_zm(drop = TRUE),
linetype = "solid", size = .4, alpha = .8, color = "tan4"
) +
geom_sf(data = units_subset, aes(fill = as.factor(burn_status_simple)), show.legend = FALSE, color = "white", alpha = 0.75) +
geom_sf(data = spatial.boundaries.parks, fill = NA, size = 0.7) +
scale_fill_manual(values = cols) +
scale_y_continuous(breaks = seq(40, 50, by = 0.02)) +
scale_x_continuous(breaks = seq(-95, -90, by = 0.02)) +
theme_minimal() +
theme(
axis.title = element_blank(),
plot.title = element_text(size = 11, face = "bold", color = "orangered3"),
plot.subtitle = element_text(size = 10, face = "plain")
) +
map_extent_park(park) +
annotation_scale(location = "br", style = "bar", bar_cols = c("grey40", "grey40"), line_col = "grey40", text_col = "grey40", line_width = 0.5) +
labs(
title = "Prescribed burns",
subtitle = paste0(park, ": {next_state}")
) # not sure why next_state works but closest_State doesn't; closest lags the data by 1 year
# Facet by year for animation frames
plot_animation <- plot + transition_states(
year,
transition_length = 0,
state_length = 1,
wrap = TRUE
)
# Generate animation
mapGIF <- animate(
plot_animation,
nframes = length(unique(units_subset$year)) + 2,
fps = 1.5,
end_pause = 2,
renderer = gifski_renderer()
)
# Save animation
anim_save(
filename = paste0("burns_", park, ".gif"),
animation = mapGIF,
path = here("Results", "990_burn_animation")
)
}files <- list.files(path = here("Results", "990_burn_animation"), pattern = ".gif", full.names = TRUE)
files_names <- list.files(path = here("Results", "990_burn_animation"), pattern = ".gif", full.names = FALSE)
files_names <- gsub("\\.gif", "", files_names)
files_names <- gsub("burns_", "", files_names)
for (i in 1:length(files)) {
cat(paste0("#### ", files_names[i], "\n"))
cat(paste0("\n\n"))
}By Sam Safran